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We present the results of a numerical study of a model of the hydrodynamics of a sheared ne- 
matogenic fluid, taking into account the effects of order parameter stresses on the velocity profile, 
but allowing spatial variations only in the gradient direction. When parameter values are such that 
the stress from orientational distortions is comparable to the bare viscous stress, the system exhibits 
steady states with the characteristics of shear banding. In addition, nonlinearity in the coupling of 
extensional flow to orientation leads to the appearance of a new steady state in which the features 
of both spatiotemporal chaos and shear banding are present. 



Experimental observations [1-5] of complex dynamics, 
including spatiotemporal chaos, in sheared wormlike mi- 
cellar solutions have stimulated the development of sev- 
eral theoretical models. The linearly extended nature 
of wormlike micelles leads naturally to considerations of 
models in which a nematic order parameter field 

is coupled to the hydrodynamic velocity. Spatiotempo- 
ral rheochaos was demonstrated [8] in the equations of 
passively sheared nematic hydrodynamics, with spatial 
variations allowed only in the gradient direction. This 
study, however, did not find any clear signature of the 
formation of shear bands observed in experiments [3- 
A different approach j9l-[l2j. based on the Johnson- 
Segalman model [13], with an added diffusive term (DJS 
model), shows shear banding but no instability of the 
shear-banded state for spatial variations only in the gra- 
dient direction. If spatial variations in the flow [10] or 
vorticity directions are allowed, the shear-banded 

state exhibits an instability that leads to complex, pos- 
sibly chaotic dynamics. 

In this Letter, we present the results of a numerical 
study of a model [8] of the hydrodynamics of a sheared 
nematogenic fluid in which spatial variations are allowed 
only in the gradient direction but the assumption of pas- 
sive advection [8] is removed, so that the effects of order 
parameter stresses on the velocity profile are fully taken 
into account. For parameter values such that the stress 
from orientational distortions is comparable to the bare 
viscous stress, the system exhibits a new steady state 
where we see the characteristics of shear banding in ad- 
dition to the states seen earlier [H. Further, allowing 
nonlinearity in the coupling of extensional flow to ori- 
entation, leads to the appearance of new steady states. 
Among these new attractors of the dynamics, the most 
significant one combines the features of both spatiotem- 
poral chaos and shear banding. Thus, going beyond the 
passive advection approximation allows the occurrence of 
banded chaotic states which were not observed in Ref. Q . 

In the simplifying limit where nonlinearities arising 
from the free-energy functional for nematic order are ig- 
nored, our model is equivalent to the DJS model studied 
in Refs. (9l-[ll|. whose results we reproduce in the corre- 



sponding parameter range. However, our fully nonlinear 
model also exhibits band formation in a different region 
of the parameter space and, in contrast to the results re- 
ported in Refs. j9l4l if . instabilities of the banded state 



are found even when spatial variations are allowed only 
in the gradient direction. These instabilities lead to a 
spatiotemporally chaotic state that also exhibits features 
of shear banding. Thus, order parameter nonlinearities 
arising from the free-energy functional play a crucial role 
in our model, leading to the appearance of new attractors 
with complex dynamics that are not present in the DJS 
model. 

We now describe the model we consider and the nu- 
merical method used in our study. The nematic order 
parameter field Q(r) in our model is a traceless, symmet- 
ric second-rank tensor whose eigenvectors and eigenval- 
ues describe respectively the directions and magnitudes 
of local anisotropy. The equilibrium behavior of Q is 
assumed to be governed by the Landau-de Gennes free- 
energy functional 
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with phenomenological parameters A, B and C governing 
the bulk free-energy difference between isotropic and ne- 
matic phases, and Fi and F2 related to the Frank elastic 
constants. In mean-field theory, the isotropic to nematic 
transition occurs when A decreases below A* = 2B^/9C. 
The equation of motion obeyed by the alignment tensor 
is 



dO 1 
^ + u.VQ = r"^G + (aoA^- 
ot 



(2) 

where the subscript ST denotes symmetrization and 
trace removal, u is the hydrodynamic velocity field, k 
= (l/2)[Vu + (Vu)^] and Q = (1/2) [Vu - (Vu)^] 
are the deformation rate and vorticity tensors, respec- 
tively. The flow geometry imposed is plane Couette with 
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velocity u = y'^x. We will refer to y and z as the 
velocity (u), gradient (V) and vorticity {uj) directions 
respectively. In Eq.(|2]), t / is a bare relaxation time 
and Q^o aiid ai are parameters related to flow alignment, 
originating in molecular shapes. Lastly, 



5F 



[AQ - V65(Q.Q)5T + CQQ : Q] 



riV2Q + r2(vv.Q)5T 



(3) 



is the molecular field conjugate to Q. The contribution 
of the alignment tensor to the deviatoric stress is given 

by 



o-oP = -aoG-ai(Q.G) ST' 



(4) 



The total deviatoric stress is cr^P plus the bare viscous 
stress. 

Under the passive advection assumption, the bare vis- 
cous stress of the system is constant and it is sufficient 
to consider cr^P alone, as was done in Ref. [8]. To incor- 
porate the full hydrodynamics in the problem, we now 
remove this assumption. We work in the Stokesian (zero 
Reynolds number) and incompressible limit, as is appro- 
priate for reasonable experimental realizations of the sys- 
tems of interest here. Thus, 



V7 ^total n 



(5) 



where cr^^^^^ is the total stress tensor in the system, and 



V.u = 0. 



(6) 



We consider spatial variation only along the ^-axis. 
Then, Eq.(|5]) with gradient terms involving derivatives 
only along the ^/-direction reduces to 



da] 



op 



dy 



(7) 



where i = x^ ja is the shear viscosity and u = yjH + 
^ix + (52Z, where 5i and S2 are ^/-dependent perturbations 
in the velocity profile. Since the fiuid is incompressible, 
and spatial variation is only along the gradient axis, per- 
turbations in Uy are zero. 

Following [3, m, time is rescaled by r/A^ and Q by 
Q/c, its magnitude at the transition temperature. Dis- 
tances are rescaled by the diffusion length constructed 
out of Fi and A^. The ratio F2/F1 is a parameter set 
to unity. We define a dimensionless viscosity parame- 
ter T] = jii/{aorQk) and consider two cases, rj = 1 and 
T] = 100. The ratio of the bare viscous stress to the 
stress from orientational distortions is determined by the 
quantity r]j where 7 is the dimensionless shear rate. 

When the right-hand side of Eq.(|4]), including nonlin- 
earities arising from the free-energy functional, Eq.([T]) 
and the ai term, but excluding the gradient terms, is 



linearized in the deviation of Q from its uniform aver- 
age value - zero, if the underlying equilibrium phase is 
isotropic - the derivatives of Q in the left-hand side of 
Eq. (|2]) are easily re-expressed as derivatives of the order- 
parameter stress. This equation, with further lineariza- 
tion of the nonlinear terms in G and a redefinition of 
parameters, then becomes equivalent to the constitutive 
equation for the viscoelastic stress in the DJS model. 
Thus, the DJS model may be thought of as a simpli- 
fied (linearized) version of the isotropic-phase limit of 
the model considered here and the effect of the additional 
nonlinearities present in our model on the behaviour ob- 
served in Refs. becomes a question of considerable 
importance. 

In our numerical study, a spatially discretized version 
of Eq.(|2]) is integrated forward in time using a fourth- 
order Runge-Kutta algorithm. For much of this study 
we work at A = and ai = 0, as in [8], so that the 
system in the absence of shear is deep in the nematic 
phase, in fact at the limit of met ast ability of the isotropic 

phase [14]. Our control parameters are = ^J^oiq re- 
lated to the tumbling coefficient in Leslie-Ericksen theory 
[7], and the dimensionless shear rate 7. The two addi- 
tional equations, Eq.([7]), that enforce the Stokes condi- 
tion, are solved simultaneously with the equation of mo- 
tion, where the matrices, k and have terms involving 
5i and 82- Starting from specified initial conditions, we 
let the equation of motion evolve and obtain the order 
parameter stress. Using cr^P, Eq.Q is then solved to 
construct the updated velocity profile, which is then fed 
back into the equation of motion. Spatial derivatives are 
approximated by symmetric finite differences defined on 
the sites of a uniform mesh of N points. Boundary con- 
ditions are so fixed that 81^82 = at the walls. Also, 
for defining derivatives at the mesh points i = 2 and 
z = - 1, we set /o /i, /at+i /at, where / is any 
variable of interest. 

We find that for 77 = 100, this model reproduces the 
various steady states ("phases") seen in [8] with only 
small shifts in the phase boundaries. This is not sur- 
prising: for large values of 77, the velocity corrections are 
small as can be seen from Eq. ([7j) • For 77 = 1 , where order 
parameter stresses are comparable to those of the sol- 
vent, and ai — 0, we find a new phase, as can be seen 
in Figs. 1 and 2. In this new phase, the steady state 
is a high-stress band with spatiotemporal or only tem- 
poral periodicity bounded by low-stress bands showing 
temporal periodicity (see Fig.2(e)). This new phase is a 
well-defined, banded periodic state between the irregu- 
lar chaotic phase and the flow- aligned fixed point. The 
position of the high-stress band thus formed depends on 
the boundary conditions. The above case is seen when 
the order parameter tensor at the two ends of the system 
is aligned in the shear plane, irrespective of orientation. 
If, however, the tensor at the two ends is aligned paral- 
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lei to the walls of the Couette cell, i.e. along then a 
low-stress band with temporal periodicity is formed be- 
tween two high-stress bands with either spatiotemporal 
or temporal periodicity near the walls. Space-time plots 
of the shear stress in the phases found for 7^ = 1, ai = 
are shown in Fig. 2. 
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FIG. 1: Phase diagram for 77 = 1.0 and ai = 0. Phase bound- 
aries between regular periodic, chaotic, banded periodic, and 
aligned (fixed point) attractors in the 7 — Afc plane are shown. 



the parameter-space region where a coexistence of high 
and low shear stress bands is first seen, the points of in- 
terfaces between the bands show complex oscillations in 
time. But as one moves deeper into this banded phase 
with increasing A/e, the complex oscillations die down to 
simple periodicity and the order parameter stress exhibits 
only regular periodic character in the banded phase [Fig. 
3(i)]. With further increase in A/c, for a constant shear 
rate, the interfaces between the high-stress band and the 
adjoining lower stress bands merge and we see only tem- 
porally periodic oscillations, with a phase difference be- 
tween the temporal oscillations at the merging point of 
the interface [Fig. 3(j)]. As A^ is increased further with 
7 kept constant, chaotic oscillations build up from this 
temporally periodic state. Obtaining a complete phase 
diagram for ai ^ is difficult, mainly because we find 
multiple locally stable attractors (and consequent depen- 
dence of the steady state on initial conditions) in some 
regions of the 7 — A^ plane. 



Periodic-->Chaotic-->Banded Periodic-->Aligned 




300 




300 


250 




250 


200 




200 


150 




150 


100 




100 


50 




50 




a) 100 200 300 




/=3.5,?t =1.1(a),1.115(b),1.135(c),1.16(d),1.36(e),1.385(f) 



FIG. 2: Space-time plots of the shear stress in the phases 
observed [(a): regular periodic, (b)-(d): chaotic, (e): banded 
periodic, (f): aligned] for a fixed shear rate, 7 = 3.5 and 
varying A/c, with 77 = 1, ai = 0, A = 0. 

For non-zero values of ai and 77 = 1, we obtain not 
only spatiotemporally periodic but also spatiotemporally 
chaotic attractors, where a high-stress band coexists with 
low-stress regions, as seen in Figs. 3(b) and 3(c) for 
ai = 0.3. The position of the high-stress band is not 
fixed, and the steady state shows such a band of consid- 
erable width along the ?/-axis (gradient direction), me- 
andering between the two walls of the Couette cell as a 
function of time. This spatiotemporally chaotic state (see 
below) with a wide high-stress band is found only for rel- 
atively low shear rates. At high shear rates (7 > 3.7), in 
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FIG. 3: Phases for 7 = 3.5 with Afc = 0.888(a), 0.945(b), 
1.025(c), 1.055(d), 1.07(e), 1.137(f), 1.19(g), 1.33(h), ry = 
1.0, ai = 0.3, A = 0. Space-time plots of the shear stress are 
shown in the regular periodic [(a)], banded chaotic [(b) and 
(c)], chaotic [(d)- (f)], banded periodic [(g)] and aligned [(h)] 
phases. Phases for 7 = 4.1 and Xk varying between 0.888 and 
1.33 are the same as those for 7 = 3.5, except for Xk between 
0.91 (i) and 1.04(j). For 7 = 4.1, the banded chaotic attractor 
of panel (b) evolves into a banded periodic attractor of panel 
(i), followed by a simple temporally periodic attractor with a 
kink [panel (j)], and the chaotic attractor of panel (e) as Xk 
is increased. 

To characterize the chaotic states found, we study their 
Lyapunov spectra (LS). For a discrete N dimensional dy- 
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namical system, the Lyapunov exponents A^, i = 1 : A^, 
arranged in decreasing order, form the LS. The number 
Nx^ of positive Lyapunov exponents and S^^, the sum 
of the positive Lyapunov exponents, are useful quantities 
that can be calculated from the LS. In particular, 
provides an upper bound and often a good estimate for 
the Kolmogorov- Sinai entropy that quantifies the mean 
rate of growth of uncertainty in a system subjected to 
small perturbations [l5|. Both these quantities scale ex- 
tensively with system size in spatiotemporally chaotic 
systems. Since our system is an extended one with a 
large number of degrees of freedom, computing the LS is 
difficult owing to the inordinately large computing time 
and memory space required. The LS of a subsystem gov- 
erned by the same equations of motion, when suitably 
rescaled, can lead to the LS of the whole system. So 
instead of analyzing the whole system, we consider sub- 
systems of comparatively small size Ns^ at space points 
j in an interval io < j ^ + J^s — 1, where io is an 
arbitrary reference point, and study the scaling of the 
relevant quantities with subsystem size Ns- For spa- 
tiotemporal chaos, it is also expected that with increasing 
subsystem size Ns^ the largest Lyapunov exponent would 
increase, asymptotically approaching its value for subsys- 
tem sizes of the order of the system size. Our analysis 
shows that the embedding dimension at certain reference 
points can be so high that the scaling with subsystem 
size can only be partially studied due to computational 
constraints. However, over the limited range of subsys- 
tem sizes that we can access, we find, depending on the 
choice of the reference point zq, clear evidence for spa- 
tiotemporal chaos of varying dimensions in the banded 
state shown in Fig. 3, panels (b) and (c). As shown in the 
top two panels of Fig. 4, both Nx^ and T^x^ exhibit an 
approximately linear increase with increasing Ns for the 
three reference points chosen. Also, the largest Lyapunov 
exponent (bottom panel of Fig. 4) shows the expected be- 
havior as a function of subsystem size. 

As mentioned above, a linearized version of our model 
is equivalent to the DJS model considered by Fielding 
and co-workers and we have reproduced their 

main results for a similar choice of parameters. In par- 
ticular, we have found that the shear-banded state found 
by them remains stable when the additional nonlineari- 
ties in our model are included. However, the values of r] 
for which we find the banded attr actors in Figs. 2 and 
3 lie outside the range in which the linearized model ex- 
hibits shear banding. Also, in contrast to the results of 
Refs. [ilill, our nonlinear model exhibits an instability 
of the banded state as Xk is increased, even if spatial 
variations are allowed only in the gradient direction (see 
Figs. 2 and 3). These results imply that the physics of 
the formation of the shear banded state, its instability, 
and the coexistence of shear banding and spatiotempo- 
ral chaos (for ai ^ 0) in our nonlinear model is very 
different from that of the DJS model. Thus, our order- 
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FIG. 4: Number of positive Lyapunov exponents (top left 
panel), sum of the positive Lyapunov exponents, (top right 
panel), and the largest Lyapunov exponent (bottom panel) 
as functions of subsystem size Ns for rj = 1.0, ai = 0.3,7 — 
3.5, Xk — 0.95. Embedding dimension for the time-series at 
each space point is different for different iq. The largest 
embedding dimension, 22, is for iq — 245, while those for 
io = 100, 400 are 18 and 20 respectively. 



parameter based model, which includes the DJS model 
as a special case, exhibits additional complex collective 
dynamics arising from nonlinearities in the free energy 
that describes the equilibrium behavior of the order pa- 
rameter. Development of similar models for rheological 
chaos in other complex ffuids would be most interesting. 
It would also be worthwhile to explore the consequences 
of allowing spatial variations in the ffow and vorticity 
directions in our model. 

We are grateful to Sriram Ramaswamy for many help- 
ful discussions. 
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